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The Kalman Filter 1s a widely used procedure in tracking algorithms. When 
normalitv assumptions are violated. the Kalman Filter performance tends to degrade. 
In this thesis a new procedure is introduced for accommodating non-normal properties 
of measurement error distributions. The procedure is developed for the multi-observer 
situation. Simulation experiment results are presented and numerical comparisons are 


made between the Kalman Ftiter performance and that of the new procedure. 
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I. INTRODUCTION 


Tracking algorithms are widely used in modern technology. Perhaps the 
most commonly used algorithm is the Kaiman Filter (KF) (cf. [Ref. 1) and] heme 
whica is based on the roilowing assumptions: 


a} The target movement model is Known, and has normaily distributed 
tluctuations. 


b) The sensors tracking the target provide unbiased, normally distnbuted 
measurements of current target “position”. 

If the assumptions are satisiied, KF is the best algorithm possibie. in practice, 
however, the normality assumption may not be justified, and the KF performance 1s 
degraded. A number of modifications of the KF have been made by different amie 
based either on heuristics. analytic derivations. or both in order to improve AF 
performance in different situations (cf. [Ref. 3] and [Ref. 4] ). 

In this thesis a tracking procedure is developed for the following one-dimentional 


basic (and classical) model: 


; = == oC ‘ ag NG ne ’ 1= i Ge 
a} ge | iS “yo iG 4. | No v rr ] () cece 

{ } ? + = t = ) ] Ze 5 f . me, : a st im = 2 3 } = { 2 D 

eats evan. eee (a); J= 123.07 7H 
where 

: Gate iS the tareet position miter time step 7 

-  y,— 1, 18 the measurement of target “position” obtained by sensor ; al ume step 


Limon ales 

in the model the target performs a random walk with normally distributed steps with 
mean zero and known variance <~. and the /” sensor's measurement errors are 
unolased ind t-distributed with known parameters (d.,5,) where d, is the degrees of 
ireedom and «. is the scale parameter of the t-distriobutuuon. There are sensors. 

This modei violates the ciassicai KF assumptions ov allowing ‘hick-talied 
‘outhier-prone} measurement errors. This model is designed for the situauon where #1 is 
small, and it isn't possible to make use of the central limit thegremme moe 

The mathematical tractability of the mode! is complicated by the fact that there 
is no simple analytic form for an estimate of tocation Of the target smith Stidemes 


MEASUrEMeENt CIrors. 
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In Chapter II two possible ways to derive the classical KF equations are given. 
These derivations will help to motivate the procedure suggested in Chapter III. Some 
practical considerations regarding the implementation of the procedure will be given in 
Chapter III as well. 

In Chapter [V the simulation experiment implementing the procedure derived in 
Chapter III is described and numerical performance comparisons of the KF and the 
Mew procedure are presented. 

hapter V contains final remarks about the new filter and further topics for 
miture development. 

Appendix A contains tabuiated results of simulation experiment. 

Appendix B contains details on simulation experiment tmplementation. 


Appendix C contains graphical presentation of the results from Appendix A. 


Mel 


It. KALMAN PIEPER 


The basic model for the classical Kalman filter is that the target moves according 
to a random walk with normally distrrbured step sizes. !he measurement errorseaee 


also normal and statistically independent lrom targee Gosiien are. 


0 
Up ae an So, Nae = wieeO 
a) ae 0. oom C+] N(O,7-)) ar = Oe 
ve = Ome : j ~\ male =e = ? oO 
o) i p+ Oe 1 eee N(O.¥,); fHi2s.u nS 


There are three alternuttve wavs to obtain tlic "Ni cuuecrome 
I. the maximum likehhood estimation approach. 
2, the minimum varianee unbiased estimation approach. 
3. the mininuzation of the mean square of estinraon error (0 — @ ). 


Extension and synergism of the frrst two methods will lead us to the procedure 


Geschncd mr iipler ada. 


A. MAXIMUM LIKELIHOOD 
A A 
Let 8 be the estimace of 8, after time step x and C, be its variance. Since 0, is 


normully! distributed with mean § and variance C,, the conditional likelihood of 
a 


7) Bs : ae ; | eae 
Ot ee { given 0 ind Yop | pr yy i: Le 
A ? 2 
io i ee ( ) 
| ntl ? 1 2 n+l 
ne 
r 0 9 \ - “) x 7 “J , A 
Ly \ eee 1 j= 2 o deed 0 
ul 1+ 1! ae ee eter t : 
} anll| 
J 


where oe by (2) = G.. nee 


Andy = oe = Van). (a= lee. 


'in the KF case normality is preserved at each estimation step as can be easily 
verified. 
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Finding the maximum for this likelihood is equivalent to finding minimum of the 
exponent. 


Le. 


| (S aot 
—— OF 


)=a0@,,, ee) 
py i 


A 
(dropping the first subscript on the ys and using onze for a compact notation). 


> 


Expanding the squares we obtain 





ani 1 Mn yo ‘rt y, 
> Sy } 
Qa .)=07. \ —-2 Yas = (2.3) 
+4 a+ v n+l — , — y 
joO ] J=O0 yo sO y 
This simple quadratic expression is mininuzed by letting 
oe 
V 
im Ee , 
aya: n Aalst 
oo ate fee 
— ee wae ey 
A y="; Cit Si J 
(K.1) 8 = = 
eo ca) 
es ——~ + See 
ey ae Saaetew 
j=0oj Co eta 


since (K.1) is a linear equation we may immediately obtain the variance of the 


estimate. 
m se 
(ama a! 
Varn Se | 7G 
7 ' ; A ye 7-0ee ! 
2) C ae ee 
eck noel es 2 7 5 m 1 
— NS), eS 
ein! — anes 'y 
y=0 "4 f oie pao 4 


eeiavions (X.1) and (A.2) are the only ones needed to inyplement the KF. 


7] 


B. MINIMUM VARIANCE 
ven yy= and y.=y ; Jala; mt is desired to get an unbiased 
Minimum variance estimate (UMVE) of ee i 
SINCE Yow yew, are all unbiased estimates of 0+ 1», any linear combination 


ne a Gy t.. Fay (2.4) 


l 
constrained by 0 1 will also be an unbiased estimate of 04 


1 
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Since 


Var({1)=a)*vy +a 


2 Zz 
gt Ov, te (2.5) 


finding the UMVE ts equivalent to solving the following problem: 


> 5 > 
paint V = PU, 
Oa | iia 


nee (2.50 
S./. @, +4 
() ] rl 


Given a Lagrange multiplier A, the Lagrangian equation may be expressed as: 


2 at > ' ' 2 “ 1 ze | 1 om 

ie Un Vo +> tL 5 oi gt aig Fae i My te hy ee Tee a al ys 2 /) 
Aer talkir urtial derivatives andesetime td iemecaue ZG ) y 
\{ter taking partrai derivatives and setting them equal to zero the resulting \sysiem 


Ole 2 intarcanaions 15. 


5 2 
hig A =!) 
2g, ~A=D 


OL eee CL sen ae 


Multiplying each equation by the corresponding I (2v,) and adding all the equations 
(OC EUNICE feslts im. 








\ = (2.9) 
m % 
Se 4 
ae v 
a | 
Substitution then leads to 
l 
Coss S 
; a (2.10) 
Ne 
oe 
t=O) 7 


So, once again the final equations are equations (K.1) and (K.2) 





m 
Su 
a j=0", 
(K.1) 9 =~H- 
a+rl m 
ee 
aoe 
Le 
ie 2) Core = minty = z 
a ih 
y=0"j 


Note that the minimum variance approach does not require normality assumptions, 


but normality assumptions make it equivalent to the maximum) likelihood approach. 


IS 


Hl. W-FILTER 


Now to return to the basre model, turget motion rs a random walk with normal 
step size, but the #1 sensors have independent Student-t aistrrbuted measurement 
enrOnS. 


The basic model: 


a) (} = ae ~ N(0.77): n= eo ce 


fees 0, Yay? Cel 
: = ' 5 . y nw tf re oar 7= o DS 
b) 54 + i 0+ | oie oF aS Ly? Oe Ly (a6); J ee clases t UME Ooh 


[t is assumed that our procedure 1s producing unbiased approximately normaily 

AN 
distributed estimates. This means that after time step 4, 9) is approximately normally 
distributed with mean 0, and vartance C,. This is the same assumption as used by 


“A 
[Ref. 3] and [Ref. 4]. Based on dirs assumption we need to construct O +} and Cm 


A. THEBESTO.., 


i. The Criterion 
A 


Based on the assumption above, the conditional likelihood of 8°. , given 9, 





and yp pee te yoy Has the form: 
Aa 
19 -9 . 
a n 
ae 
. ee ae eatee a cid ) 
a yo lee 
a Se ee age . | d.+1 — 
/> } 
> 
ae rs - 9 
i ( 22+ hy a l | 
| g a 
J / 


where e(d@.) 1s Some constant depending on the degrees-oi-ireegom parameter a: of the 
a sensor. Although it 1s possible to look for the maximum Itkehlood point, taking 
that approach has the following two interrelated potential problems. first there may 
be multiple locul maxima present. In an extreme case there may be as muny as 7A Fl 
local maxima. (See [Ref. 4] for local maxima conditions analysis when w1= 1.) Second 
im the muluple maxima situation, the global maxima niay be tall and iskimieee 
having much Irkehhood to support it. One would not like to commit all behef on the 


basis of such evidence. 
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To deal with these problems use the following (smoothing) approach will be 

. “A e . eye . ° 
used: instead of choosing 84, to ize the probability of having 8,4, in an 
infinitesimal probability element d@ around Cena, PteeeStiTULOr oon will be chosen 


foes tO maximize the probability of having 0 tl) =8oan interval 


erie al 
a N / 

[G, fe W004 ' +w], the lucer interval having a fmite (usually non-zero) predetermimed 
emer. im Other words the concern is not having an estinnite exactly “on the spor” 
but rather having it within distance w from true value of 0). , . 


Paeiiatiematicammeriissiie: approaci) wil be to find the solution to the 


equauon: 
eae A 
= ) 
maxi(@__,) = max Ei See ene dc 
Miv.1.a) ie 
OP es 


[In the equation (/V.1.a) w 1s serving as a tuning parameter and has the 
following effect on the solution: 
- When w=0 solving the equation (JV.1.a) is equivalent to finding maximum 
likelthood point. 
- When w>Q9O and finite, equation (fF.l.a) tends to dJown-weight skinny peaks 
more than fat ones. 


- When w is large enough. equation (H7.1.2) will have a unique~ solution. 
Miowever it 1S venenilly neither required nor optimal to use very large 
values for w. 


- When w approaches miinitv, any valre of Fy , will sausfy the equation. 

The idea of using a non-zero w has additional appeal since. m some pracucal 
Siemations, Occurrence of a Jarge tracking error max degrade sensor periormance for 
consecutive measurements. (Such sensor dependence on tracking performance is not 
reflected in the model.) 

2. The Solution Tecimique 
Finding a solution for equation (fV.1.a) directly involves an exhaustive search 
With numerical mtegration; not a very exciting prospect at first sight. [lowever it has 


the advantage ot algorithmic simplicity and guaranteed global maxima. 





Wess ; ; ~ . . o . 
“Since the hkelihood function is positive and asymptotically approaches zero 
when 04 , approaches = % this statement can be easily proved. 
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a. Analyticfiterative approach 
Alternatively, we may try to find an analytic solution to the equation 


(F¥.1.a) by taking the derivative of (8) 4.) with respect to 0,4, and setting it equal 








tO ZErO. 
A 
8 +w—-9 )- 
tat a) 
n+ { oe = 
Bie gi - 
rs | d.nxl 
21 ra j=l } 
a 
fv _— =—_ i. 4 2 
et ele iv nae) L 
: a. qd 
(3.2) 
ae 
(68 —wy—9 } 
1 wel 1 
- C a ie: cid .) 
" eae j 
: 1 | 
d.+1 
Lat . J 
7a _ — ae Ps 
[ ot oa Cra ee 
i+( St —_— ) — 
{ QO. le 
J e 
After setting the derivative equal to zero and performing cancellations 
> 
(8 ye ia 
1 a+ Z 
~ Cire i | 
2 t x { | = 
: i i +1 
-=1 j 
J 4 > = 
| (eee oe l | ~ 
| tie | See 
. g d. | 
3 4 ras 
(3209) 
) 
‘Q ipsa : 
a1 
ae ; 
~ oe we { 
5 ” = | i 
ee fd -{ 
es “Nn ay SS 
~ » 
folly) ~ tel Lei 
{+ 7 
Gd a 
a / 
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Next, taking logs of both sides of this equation results in: 























A ed a Z 
1 94,7 4-8,) a dj +l Pe ae I 
== E> SS: LA ail ae 
2 2 ee o. a 
a CG aU j=l J ] 
(3.4) 
ae a“ 2 ae “A ’ 
Boy ee } ~ qh 2 ego me | 
i Cite pel” As J 


Multiplying through by (—1), expanding the squares, and pulling out the common 


denonunator of the logarithm gives: 








‘ a“ “A AN A oN 
8° tw +0° +20 fe ee Ue m , - : m 
: asd +{j\ln{d a° +(y -~w=-9§ y"|— \’ Ind a 
ie c ae J a ieee) atl si ] 
Alaa j=1 J/= 
(3.5) 
As 9 A “A &W A A 
Btw +8228 wt28 w-28 8m , , 
- : PS ie yee 9. l= Sind 6 
2 —_— J ied ay a — ] 
C tt j=) ysl 


Collecting terms and rearranging leads to 








49 Ww = ead) m da~ tly ~—w—A } 
atl 2 Se iy met arl warn 
i Sis oc ie *illn mn : (3.6) 
af om 2 epee —- ; —_ i —_ { —ue § a 
ee * “ a my n'y ly 2 Cea. 
And finally, 
3 a) oN me. 
‘\ A os ae ao Ay: eno eee aay 
(W 1.5) 9,28 + \ (d -la , - 
an e Ue ae 2 
= Oat — —~ 
ae jars On a oa 


Equation (f7’.1.5) is the first order condition for a local maximum for 

equation (}}’.1.a). But perhaps a more natural way to check the first order conditions 
. A . sys e . ~ . . A 

at pomt@ 4, is to see i the likelihood function takes equal values at points 8 4) — 


and Pe, + yy, 
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Me. h 


When w approaches zero, equation (IV.1.6) approaches the form: 


(d, +1)(y 
4 = 6 + (C, yy 


j = 


n+lj ees 


(3.7) 


n+l 
are avl 


“A 
ado + -98 ») 
(The nght-most term is the lirst derivative wilterespescumre oe of the logarithmic 
fuliction. ) 

This is exactly the first-order condition for a maximum likelihood parameter value. 
Since Be, appears in the logarnthinie term of (/1.1.8), cquation (FP gig 
is difficult to solve directly. One solutron procedure is to use (+1.1.0) in a rectigaae 


N 
manner to obtain 8 


4, the empirical experience with this approach midicatessiaam 
ie | 
unless wv is relatively large, direct application of equation (17.1.5) does not converse sam 
MISsteCases, 
b. Integral evatuation approach 
Instead of using a sophisticated iterative procedure, the simpic exhaustive 


scarcn, procedure described *bulow will beaused: 


Detine 
0 (0 
min nD yt te ane Bae (3.8) 
vee ce maxtO w+ tai reed ie 
ra . . ° . ~ ”~ . 
The interval T={0. 0. includes all reasonable candidate points for 9 .,. It is 


possible to evaluate equation (fF.1.6} at a finite number of points (say 20) trom the 
interval T in accordance with the resolution reyutred, and taxe the one with Smilies 


iN 
olation of equation (47.1.0) as the estimate 8. 


The approach of pickmg the point with the mmimum violation of equation 

(11.1.6) does not guarantee even local maxima. flowever in some of the simulavions 

that were verlormed it proved satisfactory unless there were many Cases when 

kehhood function had large “flat spots” and sw was small. This generally happens 
“11 the measurement variances are farge and/or the number Of Observers Faigmian es 

‘\ procedure based on evaluation of first or tigher order Gervais. 

wil be computationally costly and in general will be more diflicult to solve because of 


local maxima and “flat spots” m the likelihood function. 





‘The “Jumps” between iterations may be very large and the estimate may reach 
tomaly wnireasonabie values. 


Therefore, the safest procedure for finding the solution of equation 
(}V.1.a) is exhaustive search over the interval T with numerical integral evaluation. 
This is the procedure used in the simulation studies reported in the next Chapter. In 
the simulation, search and integral evaluation were perfornied in a single DO-LOOP 
(See Appendix B for detalls on implementation). The sinulation was not written with 
computational efficiency in mind, nevertheless the computauonal burden did not reach 


moaeunacceptable level. 


B. VARIANCE OF 0 

[n order to keep the filter going once aos has been obtained, it 1s necessary to 
compute its variance C4, . Since equation (JV. !.a) is implicit and equation (/V.1.5) 
is not linear, there is no natural analogy to equation (A.2). 


[. Linearization approach 
. | Ras | 
One possible approach to calculate the variance or 0) , 1S {0 approximate 


27 ale 

equation (/V.1.6) with a linear equation obtained by a first order Tavlor expansion of 

Pde tOvarithimic terms around a (=3',+,, for cach jy. Cancelling and rearranging 
Sci u 


y 


terms results in the following approximation to equation (JF. 1.0): 


my 
AN A ; 
(WW. 1.c) Cue eC Oe eG ec) > k- 
n+l n n — poate 
j=) 
where 
d ri 
2 = —— (3.9) 
F doo bw” 
ens 
ee ee pia) Se (3.10) 
ead 
fens leads to: 
» me 
p= heeCGm iC st)’ kv?) (au) 
n+l n ome Pa) 
j=l 


where v= Vary + Lj 


Empirical evidence indicates that in practice this approach performs very 
badly, giving values for C, which are far too small. The C,, depends only on m and 
converges to a value for large 7. 

2. Minimization approach 
a. The Idea 

[laving rejected the above linearization approach. a different ones 
introduced bv the following motivation. 

Phe best estimarey oo , 1S SOME *cOMbinaliones! Boye Lp 
In the case of the KF it is possible to pick a linear conibination having munimal 
Variince and remain consistent with maximum likelihood ctiteria. In Ue) Wea 
this consistency no longer holds, since generally the estimate obtained by finding the 
solution to the equation (/V.i.a) wil differ from the nunimum ‘variance estumate. But in 
this cause it ts stull possible find a minimum variance lincar combination constrained 


s 


to yicld the “best” (already known) value for 0) ve 


This tneans that m order to Obtiin adm Cshiinatcue | eeieiance Cus we have 





1 
to solve the following constrained mininuzation probjenn: 
The svinbolic formulation would be 
. * 4 i 4 t Bs ~ 
pened Nay re te a 3.12) 
| . ( Dited 
0 i 
0.7 +} + ota sy = 
D4. 9, Ove i a oe Un Ly ee 
ete ee ane 
0 l mn 
where 
oy 
i ee eee (3.13) 
” ny 
e | 
: aa pSlee a (Seiae 
y = Vany | J 
J arly d —2 J 
b. Formula derivation 
[oar=t. then no nummizauon is necessary and the problem has a unique 
solution. 


A 
For notauonai sunpiteity the first subseript will be dropped, and let ‘oe 


i, 
tJ 


A A 
and @=0),,. LetA, and A, be Lagrange multipliers. The Lagrangian formulation 


is then: 


m m m 
205, -(Eo8)-fE et) 
P= Sav, -r| > ay -8} -, 29, 1 (3.15) 
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Taking partial derivatives with respect to i; and A, , A, to obtain a system of m+ 3 


linear equations with v7+ 3 unknowns gives: 


oe a = 
=n M1 Vg ho Q 
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tie i XS. = 0 (3.16) 


Multiplying each one of the first #1+ 1 equations by the corresponding I, Y, and adding 


all together results tn: 


mi oy es 
—-_ —_ L ; ae 
ee eee Coan 
% 
J} =9 iy p29 J 


Since \ a.= |. 
=e 


Multiplying each one of the first »1+ | equations by corresponding y,/v., the (t+ ne 


equation by ( — 2) and adding all together produces: 


y e ¥ “av 
ee ee =F 250 ars) 
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After changing signs a system of two linear equations follows: 


“A 
hi A+A,B= 20 
A, B+h,C=2 
i 2 
where 
2 
m on mr 7 m l 
AzN4.. ges’ 2S 
reat eae ore) 
aos JOURS 2> 0% 
Applying Cramer s rule gives 
A 
‘ 28C — 2B 
AC-B 
“A 
2A — 20B 
\ — 
“74 ae > 
7 AC — B* 
Substitution yields 
Ly 7A, 
C= Y=0,1,2,...,m) 
J 20) 
J 
or explicitly 
» 
mn 1 ne Ae rm Ae m y 
y (8 Som aS ee 
js aah) a | Tee eee 
p=) s0°j =~: =) 
l= 2 . 
J Ey Oe 1 a Sead = 
= ee 
of \ ers On acer it Wc saeee / 
J=9 oe = 0b e 
Note that the denonunator is equal to zero in case i) 21 7 


improbable event. 


24 


(3.19) 


(3.21) 


Which is a highly 


Since the Ilessian matrix 


2v9 0 ; 
() “yy 0) 
Q () ieee oT 
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IS positive definite ( ye 0; for all /), the solution is a global minimunn. 


The final equation for the variance estimate is: 
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Now the equations (FV. l.a) and (77.2) contain all that is necessary to implement the 


WI procedure. 
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Seeds Wr PROCEDURE 


; eek ae 
~ . — ? 3 7/ +y : aes. ’ ee x ai or as : . 
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Fix the actual search resolution R (depending on T=/0 ae oe): 


Set up the grid with resolution R over the interval (0. say, UPA + 1). 
geek = wR) 


4. Evaluate the likelihood function at each grid point. 


U2 


Cyr 


Find the eras by picking the “best” grid point solution of equation (IF. 1.a) in 
interval T using sum of 2k+1 adjacent likelihood function vaiues for integral 
Approxmnatuion. 

(See Appendix B for details on impiementation). 


6. Compute C pyeequation (1) 2). 


mel 


BD PRACTICAL CONSIDERATIONS 
‘A number of issues related to implemenung the WE procedure on a digital 
computer are listed below. These considerations were used in in simulation described in 
feemmext Chapter. 
1. Grid sefup 
ieieeckMdlicmvemsearch over the interval Oe eer may be performed 
over a pre-specified number of points to be checked or using a pre-specified search 


resolution. 


The first approach is wasteful when T is small and the second approach is 
wasteful when T is large. Thus, it would seem reasonable to specify both a desired 
resolution and an upper bound on the maximum number of points to be checked, 
which remain constant over time. 

2. Very dense observations 


[In the case when | is verv short there may be two complicatiome: 


5) ™ % = 4 /A e 4 ~ 

{} Because of resoluuon problems 9.) may take exact value of 0. or 0. 
In ths case equation (//.2) may force C,., equal to the corresponding i 
which is not right. 


2) Compuune © 


wey OV (T¥.2) with y, values very close together may  yteld 
underflow and division OV Zero: 

in order to prevent these problems the following modification is suggested 
and was used in the simulations: if T is very short (sav less than 3 times resoluuon) 
set 9, to be equal to midpoint of T and compute C_., using (4.2). 

5. “Very smail” numbers 

When #1 is large, the likelihood function takes very low values. To avoid 
underflow. simply multiply it by a large constant. The constant that was used in the 
sunulation is 1024”, 

One more modification was implemented to avoid arithmetic with very small 
numbers. if the exponent of the likelihood function first term was less*inaie ae 
the likelihood function was set to be equa! zero. The meaning of the modification 
1s that ine. ore cedtize considers anv Large se steg larger than 
~ 7 X (standard deviation) to be a practically an impossible step. (One must be Care! 
‘with such modifications if large target ‘umps are possible, 2s could occur 1 thewmgde! 


nvolved Student-t 9-jumps), 
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IV. SIMULATION RESULTS 


A. SIMULATION DESCRIPTION 

A simulation experiment was performed on an IBM 3033 and 4381 at the Naval 
Postgraduate School. The programming language used was FORTRAN-77, and 
single precision was used. Normai and Chi-square random numbers were generated 
with the aid of [MSL librarv routines. Student-t random numbers were generated 
from normal and Chi-square randoin varlabies with appropriate degrees of freedom 
(see [Ref. 5] ). Tracking error statistics were computed using the HISTGP subroutine 
from the NONIMSL library at NPS. Figures were generated using an IBM 
experimental graphics package. GRAFSTAT. 

For each simulation replication the tracking sequence was performed for 100 
Beeos , (le. 77 1—1,2,....100). [he number of replications was 1000 in ail cases. 

The KF and WF procedures (with diferent values for w) were carried out using 
the same random numbers. (With normaily distributed target steps and Student-t 
distributed measurement errors). - 

Statistics were collected on tracking errors 9 
fae — 1,25,50,75,100. 


The simulation was performed for 6 cases: 


te ae for 


no ace 


- Case | -one observer. o, = | 


po 


gee ase 2 -two observers, o, =0c,= 
- Case 3 -three cbservers. 5, =C,=c0.=! 
L -_ 


3 
- Case 4 -three observers, So ae ox= 5 


- Case 5 -five observers, oy ee OO O. i 
- Case § -five observers, C=O, — Gs os 5 


In all cases the measurement errors of the observers have the Student-t distribution 
with d=3 degrees of freedom and scaie o., and the target has normally distributed 
step sizes, with standard deviation += 1. in all cases the WF procedure was performed 
tor values w=0,1,2,3. Excluding case 1, the desired grid resolution. DR, was chosen to 
memommnin case |, the desired resolution was equal to 0.05. Case 1 is most susceptible 
to problems resulting from a very short interval T described in the previous chapter. 
In all cases the number of grid points on the interval T was limited by upper bound 


of 50. (See Appendix B for details on resolution and grid setup.) 
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After fixing a grid, the likelihood function was evaluated at each grid point over 
the interval of length T+2w. The integral was approximated by the sum of 2k+1 


values of the likelihood function with k being the largest integer not exceeding w/R. 


B. WF PERFORMANCE 


Figures C.1 through C.6 dispiav comparisons of the standard deviations* 


of 
the tracking error. 8. =e . for the XF and the WF for cases 1 throuigaues 
respectively. In ail the cases. the WF performance is considerably better than the KF 
pertormance. See Appendix A for tabulated resuits. 

Figures C.7 and C.$ display WF tracking error standard deviations with different 
values of w. The optimal wis somewhere in the vicinity of 1 or 2, but sensitivity to the 


exact value ot wis low. See Appendix A for exact numbers. 


C.. SENSITIVITY ANSEYSis 

For sensitivity anaiysis purposes, the tracking sequence was performed on each 
case with the same steps and measurements, but the WF procedure used the following 
(wrong) assumptions: 


- 6,=9.5 OF -0 = 1.5» onstedcmciaee aor) mon 5,= 3.0 (underestimating the 
measurement variance by a factor of 4) 


- d.=6 instead of a,= 3 (underestimating the measurement variance by a factor 
or 2) 
- t=9.5 instead of t=1 (underestimating the target step size variance bv a factor 
of +) 
Additional runs were overformed when target step size was a Student-t 
distributed variabie having 3 degrees of freedom and variance equal to lI. 
In ail sensitivity runs the SVP used w=2 and the desired resolution was equal to 
0.1. Also. the maximai number of pomts m the imterval)) 1 aseou: 
Figures C.9 througn C.l4 display WF performance under sensitivity analvsis 


conditions for cases | through 6, respectively. Ciearly the WF procedure is very robust 


“All the resuits presented in terms of sample standard deviations. ie. a square 
roots of the unbiased esumates of populauon variances. Means are not presented 
since, tn all cases for both filters. mean error values were very small. RMS (square 
Roots of Mean Square error) values were also computed in the simulation. In all runs 
inspected (the vast majority from ail runs) RAIS values were were slightly smaller than 
the standard deviations values (due to division by N in RMS computation and division 
by N— 1 in standard deviations computation). 


28 


with respect to measurement error parameter estimation and less robust with respect 
to major alternations in target movement model. See Appendix A for tabulated 
results. 

The strange behavior of the standard deviation plot with t-distributed target 
step sizes is explained bv the following phenomena: in very rare cases when the 
estimate 4 and actual target position 9 a+ are extremely distant (the target makes a 
huge step), WF starts to “ignore” the incoming measurements, regarding them as 
outliers and sticks to the old estimate. In such situation, C,, is increasing steadily but 
Slowly (by |1t each time step - due to the wrong target model). After some time 
(which may be as many as 30 time steps), C, becomes very large and WF stops 
disregarding observations and the estimate shifts toward the actual target position. 
(The last modification described in the previous Chapter was harmful in this 
situation. ) 

Such undesired behavior mav be treated by ad hoc adjustments to the WF 
procedure. For example, after detecting a steady increase in C, we may arbitrary set 
C, equal to a large number and recalculate the last estimate. [his sind of adjustment 
may be proper when actually implementing WF if there is concern about the normality 
of target step size. None of these adjustments were implemented in the simuiation 
experiment. it is anticipated that if an outlier-productive model of the target motion 
menemuscd, the effect noted would be automatically reduced, but to date no 
investigation has been made. 

Another comparison was made between the WF and KF procedures ov 
performing a simulation experiment with normally distributed measurement errors. 
The simuiation was performed for cases la-6a. Cases la-6a correspond to cases !-6., 
respectively, with normaliv distributed measurement errors having the same 
Varlance as t-distributed errors in the original cases. Once again the WEF used w= 2 
and the desired resolution was set to 0.1 with a maximum number of points in the 
iierval i ot 50. 

figures C.i5 and C.16 display the comparison of WF and KF procedures 


— 
bt 
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for cases la-6a. [t is evident that the advantage of KF over WF in cases la-6a is less 
than the superiority of WF over KF in cases 1-0. See Appendix <A for tabulated results. 

In order to see how much measurement error assumptions may be violated, 
several more runs were made using the Cauchy distribution (Student-t with 1 d.f.) 


Mieecencrating measurement errors (WF assumed 3 d.f. and the correct second 


Zo 


parameter 6). In general the tracking was steady. However in very extreme cases, 
when the interval T was extremely long and the old estimate 8 was not one of it’s 
endpoints (actual resolution in such cases became extremely large), WF produced 
an unreasonable estimate and did not recover from the huge error. 
To overcome the problem. two alterations were made to the simmMlapiam 

implementation: 

- The Maximum number or points to check was set to 500. 

- In the case when aif grid points on the interval T corresponded to the integral 


having a value equal to 0, the new estimate was set to be caual to J , the old 


estimate. The original implementation produced O nin for an estimate in such 

cases. (See Appendix B for a better approach bv altering end setup.) 
After making these two alterations. the runs were repeated. The tracking was steady 
with no special problems. The results are tabulated in Appendix A and may be 
considered yer, eood: 

A selected number of cases was chosen to perform a limited yanaomi 
demonstration. Simulation rums were made with 5 different pairs of seeds (on¢ for 
target steps and one for measurement errors) for 1,2.3 and 5 identical observers 
WF using w=2. The results are displaved in Figure C.17. Generally, vanability 
between different pairs of seeds does not exceed the variability between separa 
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V. FINAL DISCUSSION 


The WF procedure has a very general nature. It may be applied directly to 
any type of measurement error distribution. The distribution of measurement error 
must have finite variance and known density (even in tabular form). The WF filter 
performance with exotic distributions, using right or wrong distributional 
assumptions. needs to be investigated. 

(To appiv the WF procedure to models with non-normai target step size there 
must be an efficient wav to evaiuate the conditional likelihood function, or else 
additional assumptions have to be made regarding distribution of the estimate os 
before and after movement of the target.} 

When applied to models with normally distributed measurement errors 
(and normally distributed target step sizes), the WF procedure (with any w) ts 
theoretically identical to the KF procedure. in practice, nowever, there might be some 
Suirerence in performance due to finite resolution. 

WF has a built-in tuning parameter: the value of w. The optimal value for w is 
definitely not zero, but the sensitivity of the procedure to the exact value of w is low. 
wine best values of w seem to be in the range (7. 27), where t is the standard deviation 
Sree target step size. 

Instead of using w as a static tuning parameter where it is kept constant, w mav 


vs 
Ls 
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meecomputed dynamically, based on the latest observations. For example, w= [4 
Constrained to w>0.5 and w<-4 and the usuai procedure for grid setup (see Appendix 
B) produced very good results for all cases, but was not superior FOr Une. FeSuics 
Obtained with constant w. (Dynamic w results ‘vere not represented in the previous 
(iapter). 

Another modification to be explored. if it can be justified bv operational reasons. 
Pesce Of 2 NOn-s*mmetric interval in equation (?V.1.a), Le. performing integration 
With lower limit 94, —, and upper limit 9, + uw, . 


Equation (}.2) may be used not onlv as part of WF procedure. but as part of 


NN 
» 66 “9 : ms . 
any “good” procedure specifying 0, ,., and looking for C,., . 


>Optimal rule for dynamic w calculation mav be discovered by further research. 
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The most noticable difference between the KF and the WF with respect to C, 
estimation is that in the KF case C, depends only on a, converging to some value as n 
increases, and in the WF case C, is more dynamic, depending on the actual 
measurement values. The dynamic behavior of C, reflects the variable amounts of 
information received at each measurement step. This phenomenon gives aneeimsiaes 
performance measure on how weil the filter is doing. Nome of the cases checked by 
simuiation required corrective actions due to abnormal C, behavior (1e. steady 
increase or large jumps), except when target step sizes were Student-t distributed. 

‘A limited inspection of C, vaiues calculated by the WF procedure indicates that 
the C, values computed are somewhat larger than the actual variances of tracking 
errors. [t may conjectured that there is room for improvement im then 
procedure. On the other hand, forcing C,, close to its true value (known from previous 
simulations) produces good but not superior results to those obtained in the regular 
Way, using equation (?V.2). This suggests that the WF procedure is fairly robust 
with respect to modest changes in the C, computation process. (A constant value for 
C may oe used if there is sufficient knowledge of the environment, i.e. the C, values 
that should be obtained are approximately known and no large jumps in C,, values are 
expected). 

Extending the WF procedure to the multidimensional case does not seem to 
pose conceptual problems. Equation (W.1.a) may be replaced by an integral over a k- 
dimensional rectangle. The analog to equation (/#/.2) may be obtained by minimuzing 
the determinant of the covariance matrix. However, the computational burden of 
soiving equation (/V.1..2) in the k-dimensional case and minimizing the determinant 
mav oe excessive. io actually implement the WE procedure Sigg 
multidimensional environment svouid require an etficient wav to solve equation (HV.1.a) 


and to minimize the determinant of the covariance matrix. 
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APPENDIX A 
TABULATED RESULTS 


This Appendix contains tables of iracking error standard deviations. The 
tables are organized bv observers configuration. 


i tabie 1: Cne observer. 


rt 


folic 2) f wo identical observers. 
Table 3: Three identical observers. 
a 6Lablie 


Table 3: Five identical observers. 
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Three observers. one of whom is 3 times less accurate than the others. 
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oe 
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Table 6: Five observers. two of whom 1s 3 times less accurate than the others. 
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The following notation is used for specifving target step size and measurement 
Emon distributions: 
® N(v)- normal with mean 0 and variance v 
® t(d,s) - Student-t with d degrees of freedom and scale parameter o=s 
Pamen) two types of observers are involved, the parameters are given for the more 
accurate ones. 
All. cases above the double line use precisely the same target steps and 
meeaorement errors: only the filters are different. 
Filter notations: 
- K Kalman filter using correct variance 
- WF w=1 W-iuter using 5 d.f., v=1 and correct variance 
- WFw=i20 W-dlter using 3 d.fi. w=1 and assuming the scale of the tis 40 
RiUwehemiine wecutnuc Dabameter G fOr each observer 
- Whw=1,2d Wefilter using w=1 and assuming that the Student-t degrees of 
isco@omecanc 2a vather than the true value of @ 
- WFw=1.4t £x,W-filter the correct measurement parameters and assuming the 
StamG@ara deviation Of tne target step size is '2t rather than the 
(Re wae OF « 
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TABLE 3 
iaiREE TOP NTI GAL OBSERVERS 


Mvmunt Measur. Pie Track Length 
Step shader | ZS) 50 75 100 












(3,1) ee | 0.663 | 0.774 | 0.748 0.778 | 
-"- WE w=0 0.587 | 0.682 | 9.660 0.667 | 
ae 9.586 | 0.673 | 0.633 1.098 | 
== 600 0.669 : [O02 
a 0.694 | i* 
/W | 0.668 | 0. ). 0.662 
| OOS | 0). ) DiGi 
| | WF | Oi Sia V0: 0.193 | 0.820 | 
| Wi | Mon omeea Die eol720 10.700 
7 0.343 | 0.890 | 0.859 | 0.836 
NCP N( |___K fe Oreo 0704 | 0.768 | 
| Py ae ie 0.849 | 0.809 | 0.826 | 0.796 | 
| | 





TABLE 4 
Dini reborn veo ONE > TIMES WORSE 
















Measur. Rete _ Track Length. 
| PE orOr i aD 50 / 100 
N(1) t(3.1) | KF ie | O.Sey0.865 | 0.87) | 0.8% 
os fie WE w=0 0.642 | 0.793 | 0.766 | 0.762 | 0.769 | 
ae - [WE w=i 0641 | 0.778 | 0.738 | 0.734 | 0-760 | 
-". | poe eee eo | OS 0 3o | 0.716 ere | 
_ - IWFw=3 | 0.696 | 0.795 | 0.782 | 0.784 | 9.802 | 
| = ot Mae One 000 | 0 P6qN 774 | 0.778 | 0773 | 
| oe ie WEw= 22d 059 0.775 Pao (G2 | 0.768 
. VFiyy=3.%2¢ | 0.825 | 0.878 | 3 los 


| 13.38) ie WE w=> LOT 1455 | 0.309% 0.305 

N(1) (1.1) | WE w=? 0.807 | 0.963 | 0.989 | 0.937 | 0.958 | 

| Ri) (NG) KE , 0.770 | 0.914 | 0.908 | 0.900 : 9.878 
| -"- WRw=l | 9.798 | 0.958 | 0.937 | ).948 ; O.914 
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TABEE® 
FIVE IDENTICAL ODS Eiger 
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kK 
step | Error | 28 0 100 
Sn re _KF | 0.565 | 0.636 | 0.616 | 0.632 | 0.603 
a. ae WE iv=0 | O276 | OSGgumO ac: | 0.534 | 0m 
e | hie WEw=i | 0.480 | 0.535 | 0.524 | 0.53250 
-". | °. | WE w=? | 0.496 |; 0.337 | 0.530 | USas5 aie 
-". | -"- |WRiv=3 | 0.547 | 0.377 | 0.567 : 0.588 | 0.588 
| 2 | ae | WEw= 2,456 | 0.508 | 0.355 | 0.546 | 05007) Dies | 
| | te | WEw= 3a | 0397 | 0/542 | 0.338 | 0.552 | OST 
— | -"« | WFiv= 2.:2¢ | 0.081 | 0.629 | 0.640 | 0°07 7 RINiaea 
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0.645 | 0.787 | 0.739 | 9.751 ] 0.700 
VE w=0 O.57h)| 07647 =) Osos 647 | 0.3999 
-- | -- oe 0.567 | 0.64] | 08633) Gereas | 0.60 1 
“- | % I[WEw=2 | 9386 1 01645 | 0.639 | 0.654 | O.611 | 
a | -"- WF w=35 oly 0.671 | 9.670 | 0.697 | 0.649 : 
-"- -"- WFEw= 2.1206 [| 0.606 . 0.669) 0.6357 | 0.677 | 0.636 | 
am a | AN eee | ).384 |, Qodd | Y.od4 | 0.664 | 0.616 
af . | WFw=2tat | 0.760 | 0.745 | 0.761 | 0.793 | 0.749 
xf 802 | | | Z 
NG) NG) <F 0.694 | o.771 | 0.7901 0782 T ore 
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APPENDIX B 
DETAILS ON IMPLEMENTATION 


[In this Appendix the implementation of the simulation experiment is discussed. 
(There are many different ways to implement the WF procedure.) [n particular the 
implementation of the grid setup and search for maximal integral over interval with 


length 2w will be explained. 


ie) 6©GRID SETUP 
Two parameters (excluding w) control the grid setup: 
De DOR- the desired resolution (usually DR=0.1) 
2) NMAX - the maximal number of points to search (usually NMAX = 50) 


The first step is to find interval T=(0 ], where 


RO as 


N\ 
1 ae = min Oy aaees Aa ; 


min fies ree? 


9 6 
7 max art ei Voce) 


After finding T, the desired number of search points, n is computed: n=(T/DR)+1. 
Ifn> NMAX, then n 1s set to be equal to NMAX. 
iipene next step, the actual resolution R is computed: R=T)n. 
After xing R, the number of points on interval sv, Kk. is computed: k=(w/R]. (If w is 
computed dvnamicallv its value may change). 
Now the grid is fixed to have total number of n+2k+ 1 equally (R) spaced points 
foes i | Ponts On interval | itself (one on each end). 
The likelinood function values are evaluated for each grid point and stored in an 
array. ([n the same order that the grid points appear on the real number line). 

A more economical way to set up the grid «vould be as follows: 
Meme alate the likelihood function at the endpoints of interval T. If the 
hkelihood is positive at both ends. proceed as described above. [f the likelihood is 
essentially equal 0 on one or two endpoints, then eliminate the corresponding 
Sesemaiion, recompute interval {[, and repeat the procedure. (This was not 


implemented in the simulation experiment.) 


oy 


2. SEARCH FOR MAXIMAL INTEGRAL 

Once all likelihood values are stored in the array, the search for the maximal 
integral was implemented in the following way. 

Initiation step: compute sum of first 2k +1 values from the array. This sum is 


the integral approximation corresponding to 9 nin - This is the first candidate for 


ae i 

search loop: advance on the array each time adding a new element to the sum 
and subtracting the oldest one. If the new sumis larger than the largest among the 
oid ones, take the corresponding grid point as the best candidate so far and record the 
largest sum. 
After n+1 repetitions the procedure gives the best candidate among the’ = jegeme 
points on interval T (and the corresponding integral value). 

This tmplementation approach was chosen mainly because of its aigorithmic 
simplicity. Adding and subtracting REAL numbers repetitively might introduce 
some round-otf error. This is not a real problem in present case because: 

* 530 additions and subtractions is not a very large number. 
e The interest is in the grid point. not the corresponding integral value. 


* In close cases when round-off error ts reversing, the test cant be sure anyway 
because of finite resolution. 
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APPENDIX C 
COLLECTION OF FIGURES 


In this Appendix seventeen Figures presenting the WF performance (compared 


with the KF performance) are collected. 


® 


Figures C.1 through C.6 present WF vs KF comparison for cases | through 6 
Respectively, Wr performance is given for values of w=0,1.2.5. The cases 
correspond to following observer configurations: 


-  QOne observer 

- Two identical observers 

- Three identical observers 

- Three observers, one of whom 5 times less accurate than the others 

- Five identical observers 

- Five observers, two of whom 3 times less accurate than the others 
Figures C.7 and C.8 present WF performance with with different values of w 
mee ictire ©./ = cases 1.2.5 

seerigure C.o - cases 4.5.6 


Figures C.9 through C.14 present sensitivity of WF performance with respect to 
mode! violations. Figures C.9 through C.!4 correspond to case | through case 6 
respectively. 


Figures C.15 and C.16 present comparison of WF and KF performances for 
normaily distributed measurement errors. Comparison presented for cases la-6a. 
Cases la-6a correspond to cases 1-6. out have normally distributed 
measurement errors with the same variance as Student-t distributed errors in 
riginal cases. 

mee iaiclire Cala soresents cases 1a-3a 

- Figure C.16 presents cases 4a-6a 

Figure C.17 displavs variability of the simulation using different seeds. Five 
@amsor seeds are used. Variability is presented for cases 1,2,3.5. 
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Figure C.4 WF vs KF comparison for Case 4 - Three Different observers, 
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Figure C.5 WF vs KF comparison for Case 3 - Five Identical observers. 
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Figure C.6 WF vs KF comparison for Case 6 - Five different observers. 
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Figure C.9 WF Sensitivity Case | - One observer. 
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Figure C.10 WF Sensitivity Case 2 - Two [dentical observers. 
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Figure C.l2 WF Sensitivity Case 4 - Three Different observers. 
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Figure C.13 WE Sensitivity Case 3 - Five Tdenticinebseniess: 
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figure C.13 > WF vs KI comparison for Case la-3a - Normal measurement errors. 
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Figure C.16 WY vs KF comparison for Case Ja-6a - Normal measurement errors. 
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Figure C.l7 Variability of resuli: 
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